#clean up
rm(list=ls())

#libraries
library(apsrtable)
library(lmtest)
library(foreign)
library(survey)
library(xtable)
library(tseries)

#set working directory
#setwd("/Volumes/MONOGAN/DISSERTATION/partisanship/data/")

#load data
tx<-read.csv('tx9098.csv')
us<-read.csv('usMacropartisanship.csv')
ca.total<-read.csv('caTotalMacropartisanshipWeight.csv')
ca.hisp<-read.csv('caHispMacropartisanshipWeight.csv')
ca.white<-read.csv('caWhiteMacropartisanshipWeight.csv')
ca0<-read.dta('ca1994exit.dta')
tx0<-read.dta('texas1994exit.dta')

##### Exit Poll Results #####
##California##
#clean data and survey set
ca0$weight<-ca0$weight/1000
ca0$governor[ca0$governor==8]<-NA
ca0$race[ca0$race==5]<-4
ca1<-svydesign(id=~1,weights=~weight, data=ca0)

#create and print tables
#VOTE CHOICE: 1=brown, 2=wilson,8=other (NA)
#RACE: 1=white, 2=black, 3=hispanic, 4=asian, 5=other (4 & 5 merged)
ca.94 <- svytable(~governor+race, ca1); ca.94
round(prop.table(ca.94,2)*100)
summary(ca.94, statistic="Chisq")

##Texas##
#clean data and survey set
tx0$V15[tx0$V15=="OTHER: WHO?"]<-NA
tx0$V11[tx0$V11=="ASIAN"]<-"OTHER"
tx0$V15<-factor(tx0$V15)
tx0$V11<-factor(tx0$V11)
tx1<-svydesign(id=~1,weights=~V51,data=tx0)

#create and print tables
table(tx$V15)
tx.94<-svytable(~V15+V11, tx1); tx.94
round(prop.table(tx.94,2)*100)
summary(tx.94, statistic="Chisq")

##### Differences-in-Differences #####
##Hispanic Voters##
#subset and code CA
ca.hisp.sub<-subset(ca.hisp,select=c(caPartisanship.hisp.q,year.hisp,quarter.hisp),subset=year.hisp>=1990 & year.hisp<=1998)[,c(2,3,1)]
colnames(ca.hisp.sub)<-c('year','quarter','hispMacro')
ca.hisp.sub$post<-as.numeric(ca.hisp.sub$year>1995|(ca.hisp.sub$year==1995 & ca.hisp.sub$quarter>=3))
ca.hisp.sub$ca<-1

#subset and code TX
tx.sub<-subset(tx,select=c(year,quarter,txHispMacro))
colnames(tx.sub)[3]<-'hispMacro'
tx.sub$post<-as.numeric(tx.sub$year>1995|(tx.sub$year==1995 & tx.sub$quarter>=3))
tx.sub$ca<-0

#combine two subsets
hisp.sub<-rbind(ca.hisp.sub,tx.sub)

#fit model and generate predictions
hisp.did<-lm(hispMacro~ca*post,data=hisp.sub,subset=year>=1991); summary(hisp.did)
ca.pre<-predict(hisp.did,newdata=hisp.sub[1,c(2,4,5)],interval='confidence',level=.90)
ca.post<-predict(hisp.did,newdata=hisp.sub[23,c(2,4,5)],interval='confidence',level=.90)
tx.pre<-predict(hisp.did,newdata=hisp.sub[37,c(2,4,5)],interval='confidence',level=.90)
tx.post<-predict(hisp.did,newdata=hisp.sub[59,c(2,4,5)],interval='confidence',level=.90)
hisp.out<-rbind(ca.pre,ca.post,tx.pre,tx.post)
row.names(hisp.out)<-c("CA, 90-94", "CA, 95-98", "TX, 90-94", "TX, 95-98")

#draw figure
#png('hispanicDiD.png')
#postscript('hispanicDiD.eps', paper='special', horizontal=F, onefile=F, width=5,height=5,pointsize=10)
barplot(hisp.out[,1],ylim=c(0,80),ylab="Percent Democratic",main="Hispanic Macropartisanship")
arrows(x0=c(.7,1.9,3.1,4.35),y0=hisp.out[,2],y1=hisp.out[,3],col='red',angle=90,code=3)
box()
#dev.off()

##White Voters##
#subset and code CA
ca.white.sub<-subset(ca.white,select=c(caPartisanship.white.q,year.white,quarter.white),subset=year.white>=1990 & year.white<=1998)[,c(2,3,1)]
colnames(ca.white.sub)<-c('year','quarter','whiteMacro')
ca.white.sub$post<-as.numeric(ca.white.sub$year>1995|(ca.white.sub$year==1995 & ca.white.sub$quarter>=3))
ca.white.sub$ca<-1

#subset and code TX
tx.sub<-subset(tx,select=c(year,quarter,txWhiteMacro))
colnames(tx.sub)[3]<-'whiteMacro'
tx.sub$post<-as.numeric(tx.sub$year>1995|(tx.sub$year==1995 & tx.sub$quarter>=3))
tx.sub$ca<-0

#combine two subsets
white.sub<-rbind(ca.white.sub,tx.sub)

#fit model and generate predictions
white.did<-lm(whiteMacro~ca*post,data=white.sub,subset=year>=1991); summary(white.did)
ca.pre<-predict(white.did,newdata=white.sub[1,c(2,4,5)],interval='confidence',level=.90)
ca.post<-predict(white.did,newdata=white.sub[23,c(2,4,5)],interval='confidence',level=.90)
tx.pre<-predict(white.did,newdata=white.sub[37,c(2,4,5)],interval='confidence',level=.90)
tx.post<-predict(white.did,newdata=white.sub[59,c(2,4,5)],interval='confidence',level=.90)
white.out<-rbind(ca.pre,ca.post,tx.pre,tx.post)
row.names(white.out)<-c("CA, 90-94", "CA, 95-98", "TX, 90-94", "TX, 95-98")

#draw figure
#png('whiteDiD.png')
#postscript('whiteDiD.eps', paper='special', horizontal=F, onefile=F, width=5,height=5,pointsize=10)
barplot(white.out[,1],ylim=c(0,80),ylab="Percent Democratic",main="White Macropartisanship")
arrows(x0=c(.7,1.9,3.1,4.35),y0=white.out[,2],y1=white.out[,3],col='red',angle=90,code=3)
box()
#dev.off()

##All Voters##
#subset and code CA
ca.sub<-subset(ca.total,select=c(caPartisanship.q,year,quarter),subset=year>=1990 & year<=1998)[,c(2,3,1)]
colnames(ca.sub)<-c('year','quarter','Macro')
ca.sub$post<-as.numeric(ca.sub$year>1995|(ca.sub$year==1995 & ca.sub$quarter>=3))
ca.sub$ca<-1

#subset and code TX
tx.sub<-subset(tx,select=c(year,quarter,txMacro))
colnames(tx.sub)[3]<-'Macro'
tx.sub$post<-as.numeric(tx.sub$year>1995|(tx.sub$year==1995 & tx.sub$quarter>=3))
tx.sub$ca<-0

#combine two subsets
joint.sub<-rbind(ca.sub,tx.sub)

#fit model and generate predictions
joint.did<-lm(Macro~ca*post,data=joint.sub,subset=year>=1991); summary(joint.did)
ca.pre<-predict(joint.did,newdata=joint.sub[1,c(2,4,5)],interval='confidence',level=.90)
ca.post<-predict(joint.did,newdata=joint.sub[23,c(2,4,5)],interval='confidence',level=.90)
tx.pre<-predict(joint.did,newdata=joint.sub[37,c(2,4,5)],interval='confidence',level=.90)
tx.post<-predict(joint.did,newdata=joint.sub[59,c(2,4,5)],interval='confidence',level=.90)
joint.out<-rbind(ca.pre,ca.post,tx.pre,tx.post)
row.names(joint.out)<-c("CA, 90-94", "CA, 95-98", "TX, 90-94", "TX, 95-98")

#draw figure
#png('totalDiD.png')
#postscript('totalDiD.eps', paper='special', horizontal=F, onefile=F, width=5,height=5,pointsize=10)
barplot(joint.out[,1],ylim=c(0,80),ylab="Percent Democratic",main="Total Macropartisanship")
arrows(x0=c(.7,1.9,3.1,4.35),y0=joint.out[,2],y1=joint.out[,3],col='red',angle=90,code=3)
box()
#dev.off()

##Output Table##
apsrtable(hisp.did, white.did, joint.did, lev=.1)

##### MACROPARTISANSHIP MODELS #####
### Draw the Time Series###
#postscript('weightedMacropartisanship.eps', paper='special', horizontal=F, onefile=F, width=5,height=5,pointsize=10)
plot(us$macropartisanship.q,type='l',xlab='Year',ylab='Percent Democratic',lwd=3,axes=F,ylim=c(0,100))
lines(ca.total$caPartisanship.q,type='l',xlab='Year',ylab='Percent Democratic',lwd=2,axes=F,ylim=c(0,100),lty=4,col='forestgreen')
lines(x=ca.hisp$caPartisanship.hisp.q,lwd=2,lty=2,col='blue')
lines(x=ca.white$caPartisanship.white.q,lwd=2,lty=3,col='red')
axis(2)
axis(1,seq(1969,2010,by=4),at=seq(1,167,by=16))
box()
text(x=10,y=70,"US Total");arrows(x0=10,y0=68.5,x1=14,y1=65.5,length=.1)
text(x=148,y=65,"CA Total");arrows(x0=148,y0=63,x1=151,y1=58,length=.1)
text(x=148,y=88,"CA Hispanic");arrows(x0=146,y0=86,x1=140,y1=80,length=.1)
text(x=128,y=37,"CA White");arrows(x0=128,y0=39,x1=120,y1=45,length=.1)
#dev.off()

###Estimate three California models and one US model###
ca.intervention.hisp<-lm(caPartisanship.hisp.q~party.hisp+political.hisp+sentiment.hisp+l.caPartisanship.hisp+as.factor(president.hisp)+as.numeric(id.hisp>=199503),data=ca.hisp); summary(ca.intervention.hisp)
ca.intervention.white<-lm(caPartisanship.white.q~party.white+political.white+sentiment.white+l.caPartisanship.white+as.factor(president.white)+as.numeric(id.white>=199503),data=ca.white); summary(ca.intervention.white)
ca.intervention<-lm(caPartisanship.q~party+political+sentiment+l.caPartisanship+as.factor(president)+as.numeric(id>=199503),data=ca.total); summary(ca.intervention)
us.intervention<-lm(macropartisanship.q~party+political+sentiment+l.macropartisanship+as.factor(president)+as.numeric(id>=199503),data=us); summary(us.intervention)

#Breusch-Godfrey tests
bgtest(us.intervention)
bgtest(ca.intervention)
bgtest(ca.intervention.hisp)
bgtest(ca.intervention.white)

#Dickey-Fuller tests
adf.test(ca.hisp$caPartisanship.hisp.q,alternative="stationary",k=0)
adf.test(ca.white$caPartisanship.white.q,alternative="stationary",k=0)
adf.test(ca.total$caPartisanship.q,alternative="stationary",k=0)
adf.test(us$macropartisanship.q,alternative="stationary",k=0)

#differenced data for U.S.
diff.us<-cbind(us[2:167,-c(1,2,15)]-us[1:166,-c(1,2,15)],us[2:167,c(1,2,15)])
us.intervention.3<-lm(macropartisanship.q~party+political+sentiment+l.macropartisanship+as.numeric(id==197403)+as.numeric(id==197701)+as.numeric(id==198101)+as.numeric(id==198901)+as.numeric(id==199301)+as.numeric(id==200101)+as.numeric(id==200901)+as.numeric(id==199503),data=diff.us); summary(us.intervention.3)

###Plot of effects###
##Time index##
Time<-seq(-1,12,by=.1)

##Hispanic intervention effect##
m.h<-ca.intervention.hisp$coef[13]
g.h<-ca.intervention.hisp$coef[5]
pred.h<-rep(NA,131)
pred.h[1:20]<-0
pred.h[21]<-m.h
for (t in 22:131){ pred.h[t]<-m.h*(1-g.h^(Time[t]))/(1-g.h)}

#Long-term Hispanic intervention effect
m.h/(1-g.h)

##White intervention effect##
m.w<-ca.intervention.white$coef[13]
g.w<-ca.intervention.white$coef[5]
pred.w<-rep(NA,131)
pred.w[1:20]<-0
pred.w[21]<-m.w
for (t in 22:131){ pred.w[t]<-m.w*(1-g.w^(Time[t]))/(1-g.w)}

#Long-term white intervention effect
m.w/(1-g.w)

##Total intervention effect##
m.t<-ca.intervention$coef[13]
g.t<-ca.intervention$coef[5]
pred.t<-rep(NA,131)
pred.t[1:20]<-0
pred.t[21]<-m.t
for (t in 22:131){ pred.t[t]<-m.t*(1-g.t^(Time[t]))/(1-g.t)}

#Long-term total intervention effect
m.t/(1-g.t)

#Final plot
#png('intervention187.png')
#postscript('intervention187.eps', paper='special', horizontal=F, onefile=F, width=5,height=5,pointsize=10)
plot(y=pred.h,x=Time,type='l',xlab='Months',ylab='Predicted Shift in Democratic Support',lwd=2,lty=2,col='blue')
lines(y=pred.w,x=Time,lwd=2,lty=3,col='red')
lines(y=pred.t,x=Time,lwd=2)
text(x=6,y=5.75,"Hispanic");arrows(x0=6,y0=6,x1=5,y1=7,length=.1)
text(x=5,y=5,"Total");arrows(x0=5,y0=4.85,x1=4,y1=4.25,length=.1)
text(x=5,y=3.25,"White");arrows(x0=5,y0=3.05,x1=4,y1=2.35,length=.1)
#dev.off()
